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CONTRACT REPORT 


A CRITICAL EVALUATION OF VARIOUS METHODS FOR THE ANALYSIS OF 
FLOW-SOLID INTERACTION IN A NEST OF THIN CYLINDERS 
SUBJECTED TO CROSS FLOWS 


I. INTRODUCTION 


Cross flows over a cylinder or a nest of cylinders can be found in chimney 
stacks, aircrafts and missiles at high angles of attack, marine tow cables, offshore 
oil rigs and platforms, heat exchanger tube banks, and in the Space Shuttle Main 
Engine. These flows induce the structures to vibrate and the vibrating structures 
alter the flow field so that the fluid and the structural motions are strongly coupled 
together. It is imperative to understand the flow induced vibrations and to have 
predictive ability to avoid failure of structures. This report contains a critical eval- 
uation of various methods of analysis for flow induced vibrations of a nest of cylinders 
in cross flows as well as a critical evaluation of existing numerical methods for flow- 
solid interactions. 

Topics included in Section 2 are: description of discrete vibration inducing 

mechanisms for a nest of cylinders; experimental correlation equations derived through 
physical arguments or analytical modelings; recent developments in semi-numerical 
modeling of the flow-induced vibrations due to multiple number of vibration inducing 
mechanisms; and accumulated experimental data which are used as design guides for 
tube banks. 

Flow- solid interactions have been studied for the last two decades. Neverthe- 
less, application of major numerical analysis methods such as the finite difference 
method and the finite element method to these problems is still in its early development 
stage, and there exists only a limited number of publications related to these topics. 
In-depth prediction methods to be developed in the future for flow-solid interactions 
of a nest of cylinders in a cross flow, a problem which is geometrically and materially 
nonlinear in nature, would require major improvements in numerical analysis methods, 
unsteady turbulence models, and structural analysis methods including large and/or 
permanent deformations. Recent developments on these topics are discussed in 
Section 3. 

In addition to reviewing the flow- solid interaction analysis methods for tube 
banks, most recently available experimental data as well as experimental correlation 
equations are included in this report in order to provide a stand-alone guideline for 
design of tube banks. Also some of the finite element methods for flows as well as 
for structures, which may be used advantageously in the future development of a 
numerical analysis method for flow- solid interaction of tube banks, are discussed in 
some detail. 


II. CROSS FLOW INDUCED VIBRATIONS 


Vibration of a cylinder or a nest of cylinders in cross flows is induced by 
vortex shedding, fluctuating turbulent wall pressure (buffeting), fluidelastic coupling, 
jet switching, oscillating mean free stream flows, and aerodynamic noise generated in 


the flow field. In many practical cases, all of the vibration generating mechanisms 
act synergistically rather than separately. Consequently, detailed analytical modeling 
of the fully coupled flow- solid interaction is difficult. But depending on the arrange- 
ment of cylinders, flow velocities, and Reynolds numbers, a few of the vibration 
generating mechanisms usually dominate the rest so that rather simple analytical 
models can be developed. Most of the available analytical models are for these simple 
flow- solid interaction cases; yet these are also important tools in preliminary design 
analysis of tube banks. For practical design of tube banks, the experimental corre- 
lation equations, obtained by curve-fitting of experimental data, are used as design 
criteria. Each of the above vibration generating mechanisms, analytical models, 
experimental correlation equations, experimental data, and semi- numerical solutions 
are discussed below. 


2.1 Vortex Induced Vibration 


In the frontal area of a cylinder, boundary layers are formed on the top and 
bottom surfaces of the cylinder. As the fluid particles in the boundary layers travel 
toward the rear end, the linear momentum of fluid particles is dissipated in over- 
coming the adverse pressure gradient in the rear part of the cylinder. When the 
momentum is completely lost, the boundary layers separate and form two free shear 
layers which bound the wake. Since the fluid particles near the center line travel 
much slower than those near the free stream, these shear layers roll up into discrete 
vortices. The pattern of vortex shedding for a rigid cylinder changes as the Reynolds 
number, Re^, is changed [1]. Discrete laminar vortices begin to appear at approxi- 
mately Re^ = 40. The vortex flow becomes turbulent at Re^ as low as 150; the 
boundary layer flows become turbulent for Re^ >. 3 x 10®; the turbulent wake becomes 

narrower and disorganized up to Re, = 3.5 x 10®; and organized turbulent vortex 

d 6 

street is formed again for Re^ higher than 3.5 x 10 . 

For flow- solid interaction occurring when a cylinder is subjected to cross flow, 
only at low Reynolds numbers, vortex shedding is the major source of mechanical 
excitation. As soon as turbulence begins to appear in the flow field or if the free 
stream is turbulent (i.e., beginning from turbulence intensities 2 to 4 percent), the 
fluctuating turbulent wall pressure (buffeting) becomes a source of mechanical excita- 
tion which is of equal importance as the vortex shedding. As the free stream flow 
velocity is increased or decreased so that the vortex shedding frequency approaches 
the natural frequency of the cylinder, resonant vibration with large amplitude is 
produced . 


For a single row of cylinders placed perpendicular to the flow direction, vortex 
shedding is the major source of mechanical excitation only for small values of V If d, 

where V^ is the flow velocity between the cylinders and f^ is the natural frequency 

of a cylinder. As soon as Vg/f n d is increased, fluctuating turbulent wall pressure 

becomes an equally important source of excitation as the vortex shedding; and for 
Vg/f n d > 75, jet switching dominates the rest of vibration inducing mechanisms [3,4]. 


Cylinders in a closely spaced tube bank do not respond as single cylinders; 
rather, interaction with the flow field causes coupled motion in the group of cylinders. 
The vortex shedding is confined to a range of moderate to large cylinder spacing 
ratios in tube banks [2]. For S t /d < 1.5, vortex shedding would be suppressed due 
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to proximity to the cylinders [ 1] , where is the cylinder spacing in normal to the 

flow direction as shown in Figure 1. Idealized response of a cylinder in a nest of 
cylinders (tube bank) to increasing flow velocity is shown in Figure 2 [5]. It can be 
seen in Figure 2 that the resonance to fluid periodicity, i.e., resonance to fluctuating 
turbulent wall pressure or vortex shedding, and fluidelastic instability are the most 
important mechanisms that can cause failure of tube banks. 





Figure 1. Parameters used for a nest of cylinders. 


Vrms 


V 

Figure 2. Idealized response of a cylinder in a nest of cylinders 

subjected to cross flow. 
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Vibrations due to vortex shedding and buffeting are classified as forced vibra- 
tions since these two vibration inducing mechanisms do not require flow- solid inter- 
actions and can exist for flow around rigid bodies. The experimental correlation 
equations, used by tube bank designers to avoid resonance to vortex shedding, are 
usually developed on the basis of the contribution of Chen [6]. On the other hand, 
correlation equations developed on the basis of the theory formulated by Owen [7] 
provide design means to avoid resonance to fluctuating turbulent wall pressure 
(buffeting). Later, Paidoussis [8] showed that the two seemingly different conditions 
for resonance represent the same physical phenomenon and the only difference depends 
on the physical interpretation of the mechanism that induces resonance. In fact, as 
shown in Figure 2, the vibration due to vortex shedding, if it exists, and buffeting 
coexist in the low free stream velocity region. 

The correlation equation derived by Chen [ 6] , for resonance in vortex induced 
vibration , is given as : 


f d/V = 
n g 


vs 


( 2 . 1 ) 


where S 


= f d/V is the Strouhal number, 
vs vs g 

and the rest of the notations are the same as before 


f is the vortex shedding frequency, 

A table of S for various 
vs 

geometrical arrangement of cylinders can be found in Reference 6. Resonance to 
vortex shedding can be avoided by rendering the left hand side of equation (2.1) to 
be different, by 40 to 50 percent, from the Strouhal number which will cause resonance 
for a particular arrangement of cylinders. 


2.2 Fluctuating Turbulent Wall Pressure 

Turbulent flows induce vibration of structures through the fluctuating turbulent 
wall pressure acting on the surface of structures submerged in the flow field. 

A nest of cylinders subjected to cross flow exhibit a peak amplitude at a flow 
speed whose dominant frequency of fluctuating turbulent wall pressure spectrum coin- 
cides with the natural frequency of the submerged cylinders. Most of the experi- 
mental correlation equations used for predicting resonance which is caused by buffet- 
ing are based on the expression developed by Owen [7]. 

The correlation equation due to Owen [ 7] is based on the assumptions that : the 

conditions for fully turbulent flow have been developed through regions sufficiently 
deep inside the tube bank; the rates of production and dissipation of the turbulent 
kinetic energy are in equilibrium ; the length scale of the energy containing eddies is 
of the same order of magnitude as the cylinder spacing in the flow direction; the 
source of vibration is associated with the randomly fluctuating pressure forces imposed 
on the cylinders by the turbulent eddies; the flow velocity is low enough so that the 
fluidelastic vibration would not occur, yet the Reynolds number is high enough so that 
the effect of molecular viscosity is negligible; the tube spacing normal to the flow 
direction, S^/d, is appreciably greater than unity so that a cylinder preserves some 

of the features of a circular cylinder in isolation but not large enough for these 
features to dominate; and the rate of dissipation of turbulent energy and the work 
done by the mean pressure gradient across any row of cylinders are in balance as has 
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been observed in experiments. The correlation equation due to Owen [7] states that 
buffeting resonance would occur if: 


f n d/v 


e 



( 2 . 2 ) 


where V is defined to be the flow velocity at the minimum gap in a row, = f^d/V 
is a Strouhal number, f^ is the dominant frequency of the turbulent buffeting spec- 
trum, and the rest of the notations are the same as before. According to Owen [7], 
f b is given as: 


f = —a. 
d % 


(1 


1 

- d/S 


[3.05 


t J 


(1 - f > 

S t 


+ 0.28] 


where l is the length scale of the energy containing eddies, and the rest of the nota- 
tions are the same as before. It was mentioned that equations (2.1) and (2.2) repre- 
sent the same phyiscal resonance conditions at low free stream velocities [5]. Never- 
theless, there exists significant quantitative difference between these two equations if 
they are plotted on the same graph. Therefore, equations (2.1) and (2.2) can pro- 
vide only rough criteria for avoiding resonance to forced vibration in the design of 
a tube bank. 

A prediction method, which utilizes the random vibration theory [9,10] for the 
maximum amplitude in resonant forced vibration was developed by Chen and Wambsganss 
[11,12]. In this method, it is required to know the spectra and frequency composition 
of fluctuating turbulent wall pressure. For preliminary design analysis, neither the 
distribution of turbulent wall pressure as function of time nor the spectra and fre- 
quency composition are known a priori. Therefore, a set of data obtained from flat 
plate measurements have been used in References 11 and 12 to predict the maximum 
amplitude. The specific example problem in References 11 and 12 will not be discussed 
herein. General procedures for flow induced random vibration analysis can be found 
in References 9 and 10. 

Recently, methods for predicting fluctuating turbulent wall pressures [13,14] 
and more detailed experimental data [15] have begun to appear in the open literature. 
Use of such a prediction method in a finite element analysis of compliant plates has 
been reported in References 16 and 17. These will be discussed in Section 3.2. 


2.3 Fluidelastic Instability 

Tube banks subjected to increasing flow velocity begin to vibrate with large 
amplitude at a critical flow velocity with a definite frequency or a group of frequen- 
cies which are different from that of the natural vibration mode of the cylinders. 
Above the critical flow velocity, the amplitude of vibration increases as the square of 
the flow velocity, apparently without limit. In this type of vibration, called fluid- 
elastic or aeroelastic vibrations, the flow field excites the structure to vibrate and 
in turn the vibrating structure disturbs the fluctuation inducing flow field. Accord- 
ingly, this is a fully coupled flow-solid interaction phenomenon and is classified as a 
self- excited vibration. 
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Fluidelastic instability is of utmost concern in tube bank designs because it can 
cause total failure in the structural system. On the other hand, resonance to forced 
vibrations, i.e., vibrations due to vortex shedding and/or buffeting is of less con- 
cern, because it can be avoided by keeping the flow velocities away from the expected 
resonance velocity by a factor of 40 to 50 percent. 

The analytical expression for the onset of fluidelastic instability for a row of 
cylinders placed perpendicular to the flow direction was first proposed by Connors 
[4] , and is given as: 


V g /f b d = K fe (m<5/pd 2 )a 


(2.4) 


where the gap velocity V is defined as V = V[(S,/d)/{(S, /d) - 1 }] ; K f is a constant, 

or a function of cylinder spacing to be determined from experiments; a is a constant 
to be derived from theoretical analysis, yet it can be modified to fit experimental data 
more precisely; p is the density of the fluid; m is the added mass defined as m = m„ + 

2 U 

C up d / 4; m n is the mass per unit length of the cylinder; C is an experimentally 

determined inertia coefficient for displaced fluid due to structural vibration [ 1] ; it = 
3.14159; <5 is the logarithmic decrement of amplitude due to fluid and structural 
damping which is dimensionless; and the rest of the notations are the same as before. 
Equation (2.4) was derived by assuming that the energy supplied to the structural 
vibration would be equal to the energy dissipated by the damping force to maintain 
equilibrium in fluidelastic vibration. 

Following the line of Connors [4], Blevins [18] formulated an analytical model 
which is more rigorous mathematically. The underlying assumptions in his modeling 
are : the vortex shedding frequency of the cylinders is well above the natural fre- 

quency of the cylinders so that the cylinders are not excited by vortex shedding or 
buffeting; flow velocity is low enough so that jet switching is inactive; and the forces 
due to asymmetry of the flow pattern caused by the displacement of cylinders from 
equilibrium is the only source of vibration. Blevins [18] sub-divided the fluidelastic 
instability criterion into three cases depending on vibration frequencies and damping 
forces along and transverse to the flow direction. These are: a symmetric frequency 

and symmetric damping force case, an asymmetric frequency and asymmetric damping 
force case. The first case of Blevins [18] reproduced exactly the same result as that 
obtained by Connors [4]. Later, Heinecke [19] showed that the same equation, equa- 
tion (2.4), can be derived by simple dimensional analysis and energy consideration of 
a dampled, hormonically excited oscillator. From the underlying assumptions of 
Blevin's [18] and Heinecke's [19] models, it might be suspected that the Connors' 
criterion for fluidelastic instability is not conservative. In fact, a list of failure cases 
of tube banks, which satisfied the Connors' fluidelastic instability criterion, can be 
found in Reference 8. 

Numerous experimental correlation equations have been derived following the line 
of Connors' work by incorporating more experimental data into equation (2.4). In 
order to provide a stand-alone guideline for tube bank designs, different values of 
the coefficient and exponent a in equation (2.4), proposed by many different 

groups of workers in the area, are summarized in Table 1 [20]. 
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TABLE 1. VALUES OF K fe AND a [20] 


Investigators 

K fe 

a 

Remarks 

Connors [ 4] 

9.9 

0.5 

Tube row with S t /d = 1.42 

Blevins [ 18] 

2( 2tt) °* 5 

(C c ) 0 - 25 
x y 

0.5 

C and C are fluid - 
x y 

elastic- stiffness force 
coefficients 

Y. N. Chen [21] 

3 Re" 0 ' 25 

0.5 

Re = Reynolds number, 
3 = constant 

Gross [22] 

4/k 

1.0 

For square array, and k 
determined from fluid 
force 

Gorman [23] 

3.3 

0.5 

Suggested design guide- 
line 

Savkar [24] 

4.95 (S t /d) 2 

0.5 

For triangular arrays 

Connors [ 25] 

0.37 + 1.76 S t /d 

0.5 

For square array 
1.41 < S t /d < 2.12 

Pettigrew et al. [26] 

3.3 

0.5 

Suggested design guide- 
line 

Weaver and Grover [27] 

7.1 

0.21 

Rotated triangular array 
S t /d = 1.375 

Chen and Jendrzejczyk 

2.49 to 6.03 

0.2 to 1.08 

For various rectangular 
arrays and mixed array 
in water flow 

Tanaka et al. [29] 

3.0 

0.75 

For square array, 
S t /d = 2.0 


Another type of empirical correlation equation has been developed by means of 
similitude analysis [27]. This equation is given as: 


V 5 

S - r ( m ^ 
F"d “ B (_ 72 ) 

n pd 


(4 - i) 


(2.5) 


where a, 3, and y are exponents to be determined from experiments; B is a constant, 
or might be a function of other nondimensional parameters; the gap velocity, V , is 

b 

defined in the same way as in equation (2.4); and the rest of the notations are 
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the same as before. The experimentally determined coefficients by many researchers 
are summarized in Table 2 [20]. 


TABLE 2. VALUES OF COEFFICIENTS FOR EQUATION (2.5) [20] 


Investigators 

B 

a 

3 

y 

Remark (Based on 
published data) 

Paidoussis [8] 


0.5 

0.25 

H 

Includes all data 

Paidoussis [30] 

5.8 

0.4 

0.4 

■ 

Excludes some data 


All these theories and data appearing in Tables 1 and 2 have been assembled 
and presented by Chen [20] for different tube patterns, i.e. , according to different 
geometrical arrangements of cylinders. These are shown in Table 3. Blevins [31] 
pointed out, however, that any dependence of this data on tube pattern is masked by 
the data scatter and suggested to use a single least squares curve- fitted expression 
as shown in Figure 3 [31,32]. A similar single graphical representation of most of the 
experimental data that are available to date for fluidelastic instability prediction of 
triangular arrays of cylinders, has been compiled by Chen [33], and is given in 
Figure 4. 

TABLE 3. LOWER BOUND ON CRITICAL FLOW VELOCITY [20] 


Array 

0 

Parameter 

Range for <5 m 

V /f d 
g n 

Tube row 

1 

0.05 

< 

6 m 

< 

0.3 

1.35 (S t /d - 0.375)6m'° 6 



0.3 

< 

6 m 

< 

4.0 

2.30 (S t /d - 0.375) 6^- 5 



4.0 

< 

5 m 

< 

300 

6.00 (S t /d - 0.375) 6^- 5 

Square 

0 

0 

01 

0.03 

< 

6 m 

< 

0.7 

2.106 m 15 



0.7 

< 

6 m 

< 

300 

2. 356 0-50 

Rotated square 

45° 

0.1 

< 

6 m 

< 

300 

3.54 (S t /d - 0.5)6°’ 5 

Triangular 

30° 

0.1 

< 

5 m 

< 

2 

3.58 (S t /d - 0.9)6m 1 



2 

< 

6 m 

< 

300 

6.53 (S t /d - 0.9)60-5 

Rotated triangular 

60° 

0.01 

< 

6 m 

< 

1 

2- 86 S; 17 



1 

< 

6 m 

< 

300 

2.86 °-5 

m 


* 6 is defined in Figure 1. 
** 6 m = m6 /pd2 ' 
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Figure 3. Mean and lower boundary curves for fluidelastic instability [31]. 



Figure 4. Stability diagram for triangular arrays of cylinders [33]. 
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2.4 Other Vibration Inducing Mechanisms 

Other vibration -inducing mechanisms are: (1) aeroacoustic noise; (2) swirling 

flows; and (3) oscillating mean free-stream velocity. These vibration sources have 
received less attention than the others considered earlier in this section, and are 
more difficult to treat, both analytically as well as experimentally, especially for tube 
banks. Another vibration inducing source is jet switching which is a type of fluid- 
elastic vibration and can be observed most prominently for a single row of cylinders 
subjected to cross flows. Jet switching is included in this section. 


2.4.1 Aeroacoustic Noise 

Aeroacoustic sound is generated by vortex shedding and turbulence. If the 
vortex shedding frequency coincides with the acoustic frequency , then the two 
mechanisms reinforce each other. Furthermore, if these two frequencies approach 
the natural frequency of cylinders, then large amplitude vibration which can lead to 
failure of tube banks can be caused. 

A review of numerous analytical and finite difference numerical solution of 
acoustic wave propagation in both homogeneous and inhomogeneous medium can be 
found in Candel [34], The acoustic wave propagation equations considered in Refer- 
ence 34 are a direct wave propagation equation and the Helmholtz equation. The 
direct wave propagation equation is given as : 



c 2 V 2 P = 0 
s 


( 2 . 6 ) 


where p is the sound pressure, t is time, c is the speed of sound, c = c(x), x is 
s 2 

the spatial coordinates, and v is the Laplacian operator. The Helmholtz equation can 
be obtained from equation (2.6) by setting 



where p m (x) is the mode shape function, e is the natural base of logarithm, i = /-T, 

and a) is the radian frequency. A finite element analysis of Helmholtz equation for 
optimal design of acoustic chamber can be found in Bernhard [35]. 

A review of sound induced by vortex shedding from cylinders is given by 
Blevins [ 36] . The governing equation considered in the review is the Lighthill's 
equation given as : 



c 


2 


V 


2 


P 





( 2 . 8 ) 


10 


2 

where T~ = pv.v. - t- - c p6~, is the fluid stress tensor, 
delta, and the rest of the notations are the same as before. 


i] 


is the Kronecker 


Recently , Quinn and Howe [ 37] presented an analytical solution for acoustic 
wave propagation through a tube bank. They modeled the cylinders as rigid finite 
flat plates at zero angle of attack. The effect of vortex shedding was incorporated 
in their equation by inserting aeroacoustic dipole sources in the forcing function term 
of the classical velocity potential equation for acoustic wave equation. The result of 
this analysis showed qualitatively that both attenuation and sound speed increase with 
decreasing frequency in such a way that onset of resonance is delayed. 

Blevins [36] concluded that the degree to which self-sound synchronizes vortex 
shedding from various tubes is not known and that the role of acoustic damping had 
not been well understood until now. 


More recently , an experimental study to eliminate aerodynamic noise in a tube 
bank was reported by Zdravkovich and Nuttall [38], Experiments were performed 
with a tube bank of three rows of cylinders, with seven cylinders per each row. 

The array of cylinders is shown in Figure 5. It was shown that the acoustic 
resonance could be eliminated either by spacing unequally two successive rows of 
cylinders or by removing certain cylinders from the tube bank. Removal of a few 
cylinders from the first row was found to be the most effective way to eliminate 
acoustic resonance. For the second row, removal of a few cylinders was effective 
only at the pressure nodes; while removal in the third row was not effective at all. 

Other approaches to reduce aeroacoustic noise was also found by Baird [39] and 
Cohan and Dean [ 40] . They recommended to use baffles in order to avoid the first 
and higher modes of standing waves, but this technique may not be used advan- 
tageously for densely populated cylinders because of high pressure drop in the tube 
bank [40]. 



Figure 5. Possible tube arrangements to reduce 
aeroacoustic noise [ 37] . 
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2.4.2 Jet Switching 


In this phenomenon, fluid forces excite the cylinders to vibrate and the dis- 
placement of cylinders cause the jet pairs emerging from closely spaced cylinders to 
switch. Jet switching is a self-excited vibration which requires fully coupled flow- 
solid interaction. Coupling of jets in the wake region of a row of cylinders is shown 
in Figure 6. 



Figure 6. Coupling of jets in the wake region of a row of 
cylinders subjected to cross-flow [1]. 

It was claimed in Reference 3 that jet switching may not occur for V^/^d < 75, 

and it was reported in Reference 4 that jet switching was not observed for V /f d < 
30. S n 

It was observed experimentally that the flow slightly inside the tube bank is 
suppressed so that jet switching is confined to a few cylinders at the flow inlet 
region. The cylinders at the flow inlet region are subjected to the greatest fluid 
force as well as to jet switching. Consequently, these cylinders are most susceptible 
to damage due to fatigue and inelastic deformation due to large displacements. Closed 
form analytical solutions cannot take this fact into account. Also numerical analysis 
for this type of problem has not been attempted as yet. However, experimental 
correlation equations include this phenomenon implicitly, for experiments are concerned 
with failure of tube banks in an integrated sense. 
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2.4.3 Oscillating Mean Flows and Swirling Flows 

Oscillating mean flows and swirling flows excite a nest of cylinders submerged 
in the flow field to vibrate. Analytical analysis of flow induced vibrations due to 
oscillating mean flows is mostly based on the approach proposed by Morison [41], in 
which a cylinder is represented by an Euler beam and the fluid force acting on the 
cylinder is modeled as forcing function in the Euler beam equation. Most of the 
recent analyses have concentrated on modeling the fluid force due to oscillating mean 
flow, but flow -solid interaction has seldom been considered. 

It is natural to expect the swirling flows would induce vibration of cylin- 
ders submerged in the flow field. More importantly, the swirling component of the 
flow can alter the entire flow field so that the vibration mode of submerged cylinders 
may be different from those without swirl velocity. For example, it is well known 
that the flow field of jets with swirling component is quite different from those without 
swirling component [ 42] . Swirling flow through a nest of cylinders can be found in a 
swirl crossflow Stirling engine heater [43] or in the Space Shuttle Main Engine (SSME) 
Main Injector Assembly (MIA) (Fig. 7). A finite difference numerical analysis of the 
flow field inside MIA showed that there exists the swirling flow [44] (Fig. 8). As a 
remark, the vibration of cylinders has not been considered in Reference 43 or 44. 

Only the heat transfer rate was considered experimentally in Reference 43. Neither 
the details of flow field inside tube banks nor the effect of swirl flow on vibration 
of cylinders are known yet. 


2.5 Semi- Analytical Solutions Including Combined Excitation Mechanisms 

The vibration of cylinders in a tube bank is different from a collection of 
vibrating cylinders in isolation. More complicated mathematical models to account for 
coupled vibration of cylinders in tube banks were proposed by Blevins [ 18] , Chen 
[45], and Price and Paidoussis [46]. In all of these models, cylinders were repre- 
sented by Euler Beams and the interaction between vibrating cylinders was included 
in the governing equations by specifying the forcing functions to linearly depend on 
the displacements and velocities of cylinders. It was also assumed in these models 
that the amplitude of vibration is small enough so that the flow field is not disturbed 
due to vibration of cylinders. These models successfully describe the coupled vibra- 
tion mode between cylinders, but can not be used with confidence for problems with 
large amplitude vibrations which lead to mechanical wear, fretting, corrosion, and 
fatigue failure of tube banks. 

Blevins [18] considered ’’whirling" of a row of cylinders subjected to flows with 
low reduced velocities so that whirling vibration of cylinders is mainly due to lift 
force (caused by the asymmetry of flow pattern due to displacement of cylinders) and 
drag force. Other vibration inducing mechanisms were excluded from his analysis. 
The flow velocity for the onset of whirling instability could be obtained from the 
analysis. 

Chen [45] considered vibration of a row of cylinders subjected to a multiple 
number of excitation forces such as vortex shedding, fluidelastic coupling, drag and 
lift forces, and fluid inertia coupling. Apparently his model can be extended to 
vibration of tube banks with multiple rows once the fluid force coefficients for all the 
cylinders in the tube bank are provided by experiments or by other means, for 
example, computational fluid dynamics. The mathematical model due to Chen [45] is 
discussed below, for its generality and completeness. 
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Figure 7. Space Shuttle Main Engine (SSME) power head assembly. 



F: FUEL SIDE. OX: OXIDIZER SIDE 




0 = 20 ° 


Figure 8. Flow field inside the SSME Main Injector Assembly [44], 
(a) Horizontal velocity profile on the middle plane, (b) Vertical 
velocity profile at 30 deg plane from center line. 


The equation of motion for a single cylinder, which has the same flexural rigidi- 
ties in both flow direction and transverse to the flow direction, can be written as [45] 


E i h 


9 4 u. 
21 
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+ c i TT + m i 


aV 

9 1 2 


= f. 


(2.9) 


where the subscript i denotes the ith cylinder in the tube bank, l is measured along 

T 

the cylinder axis, E is the modulus of elasticity, I is the area moment of inertia, u = 
(UpUg), u^ is the displacement in the flow direction (hereafter denoted as x-coor- 

dinate direction), u 9 is the displacement in the transverse direction (hereafter denoted 

2 rp 

as y-coordinate direction), f = (f ,f ), f and f are forcing functions. The major 

^ x y x y 

difficulty in flow- solid interaction analysis lies in determining the forcing functions 
which depend on the displacement of cylinders. Chen [45] decomposed the forcing 
function as: 
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(2.13) 
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In equations (2.10) through (2.14), _f. c is the inertia force of fluid displaced by 
vibration of cylinders, M' is mass of displaced fluid, a., is added mass coefficients for 
ith cylinder due to displacement of jth cylinder; is the column vector of drag force 
in x- direction and lift force in y- direction, c^ £ is the lift coefficient, c^ is the drag 
coefficient, e- is a small number to be determined experimentally, w , w = 27rS/d, is 
the vortex shedding frequency, i|>. £ and are phase delay angles for lift and drag 
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forces, respectively; f. e is a column vector of fluidelastic forces, y.. is the fluid- 
elastic force coefficients; J. is the hydrodynamic damping forces, s,. and t.. are 

O J J 

hydrodynamic damping coefficients ; and f. is a column vector of forcing functions 


due to the rest of vibration inducing mechanisms. Equations (2.9) and (2.10) are 
quite general so that response of cylinders due to a specific vibration mechanism or 
any set of combined mechanisms can be analyzed by simplifying the load term, equa- 
tion (2.9). But the quality of analysis will depend on the physical coefficients used 
in equations (2.9) and (2.10). A set of analytical solutions, using the Eigen value 
analysis method, for natural frequencies and mode shapes of coupled vibrations, 
critical flow velocities for different excitation forces, and response to vortex shedding, 
lift, and drag forces for a row of cylinders have been presented in Reference 45. 


Prices and Paidoussis [46] considered vibration of two rows of cylinders in cross 
flow. Fluid force ciefficients used in the analysis were obtained from wind tunnel 
tests [47,48], The form of governing equation of motion used in their analysis is 
the same as that of Chen [45]. But they elaborated on modeling the effect of fluid 
retardness and time delay. The results showed that: analysis of a single flexible 

cylinder in a tube bank of rigid cylinders yielded reasonably good approximation of 

2 

critical flow velocity for tube banks with all flexible cylinders for m6 /pd less than 300, 

2 

but not for m6/pd greater than 300. Also, based on their analytical analysis, they 

2 

recommended to separate <$ from m/pd in the experimental correlation equations, as 
in equation (2.5). 


III. NUMERICAL ANALYSIS METHODS FOR 
FLOW-SOLID INTERACTIONS 


Geometrical and material nonlinearity inherent to flow- solid interaction phenomena 
precludes the possibility of obtaining closed form solutions to these problems except 
for simplified cases. It would be desirable to predict the flow- solid interactions 
numerically by solving unsteady turbulent flow equations, i.e. , unsteady Navier- 
Stokes equations with appropriate turbulence models. By obtaining the pressure dis- 
tribution and skin friction on the surface of the structure, one could then solve 
structural dynamics equations for vibration (and stability analysis) iteratively. The 
vibration inducing mechanisms such as vortex shedding, fluctuating turbulent wall 
pressure, and jet switching may be obtained as part of numerical solution. 

Significant progress in numerical modeling of flow-solid interactions, turbulence 
models, and nonlinear structural analysis methods has been made recently. Existing 
literature covers; unsteady turbulent flows around rigid bodies and/or rigid con- 
tainers; fully coupled laminar flow- solid interactions; and nonlinear stability, vibration, 
and elasto- viscoplastic analysis of structures subjected to known applied loads. For 
turbulent flow- solid interaction cases, the applied loads due to fluid flow have been 
obtained by solving the turbulent flow equations separately or by incorporating 
analytical or experimental correlation equations for fluctuating wall pressures induced 
by turbulence. Fully coupled numerical analysis of flow- solid interaction in cross 
flows, either laminar or turbulent, for a nest of cylinders has not been made yet. 

Some of the numerical analysis methods for laminar flow-solid interactions, tur- 
bulent flow- solid interactions, mathematical models of turbulence, and structural 
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analysis methods which may be relevant to developing numerical methods for a class 
of flow- solid interaction of a nest of cylinders subjected to cross flows are discussed 
below . 


3.1 Laminar Flow-Solid Interactions 


Most of the recent investigations using numerical analysis for flow- solid interac- 
tions have been intended for safety analysis of nuclear reactors exposed to transient 
fluid loading. In these problems, displacements excursed by fluid particles and 
surrounding or submerged structures are of the same order of magnitude, and the 
time period of interest is very short (a few milliseconds for some cases). The effect 
of turbulent viscosity has not been included in these models. Researchers in this 
area developed a special mathematical model and numerical analysis methods which are 
suitable for this special class of flow- solid interaction problems. The governing equa- 
tions are described in Lagrangian-Eulerian coordinate system which was originally 
proposed by Noh [49] and further developed by Hirt [50] to be used for gill flow 
velocities. The advantages of the Lagrangian-Eulerian description are: excessive 

distortion of computational grids can be avoided; and grid points can be displaced 
independently of the fluid motion so that implementation of fluid- structure interface 
condition can be simplified. 


Let v,(j=l,2,3) be the velocity vector of a fluid particle, v.° be the grid 
1 D 

velocity vector, and v. be the difference velocity vector defined as: 
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where j denotes the spatial coordinate directions. Then the governing equations for 
the fluid flow in the Lagrangian-Eulerian coordinate system are obtained to be: 
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where equations (3.2) through (3.5) represent conservation of mass, linear momentum, 
energy , and constitutive equations respectively , p is the density , Oy is the fluid 
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stress tensor, b. is the component of the body force along the ith coordinate direc- 
tion, e is the internal energy, p is the molecular viscosity, and p is the pressure. 

Finite element representation of equations (3.2) through (3.5) are all based on 
the same method of weighted residuals, but detailed solution techniques used by 
different groups of workers, i.e. , Donea et al. [51], Belytschko et al. [52,53], and 
Liu [ 54] , are somewhat different . 

In the finite element method, each variable is interpolated using shape functions 


nodal values of these variables. 

These are: 


Vj(x,t) = v. k (t) <f> k (x) 

k = !, N k 

(3.6) 

v. G (x,t) = v. k G (t) <j> k (x) 

k = !, N k 

(3.7) 

p(x,t) = P m (t) f> m (x) , 

m = !, N m 

(3.8) 

pe(x,t) = (pe) m $ m (x) , 

m = !, N m 

(3.9) 

p (x>t) = P g (t) x s (x) , 

s = !, N s 

(3.10) 


where the shape functions <j>, (x) , $ (x), and x ( x ) are functions of spatial coor- 

k m s 

dinates; nodal values of these variables are functions of time; and N, , N , and N 

k’ m’ s 

denotes number of these variables for each element. In general, the order of inter- 
polating polynomials and inter-element continuity requirements are different depending 
on the differential equations to be solved and the numerical method employed. The 
finite element system of equations are obtained by multiplying the differential equa- 
tions with appropriate test functions, substituting the continuous variables with 
interpolated expressions, integrating over an element, and assembling the resulting 
element system of equations into a global system of equations. 

Donea et al. [51] used linear shape functions for the velocities, and constant 
elements for density and internal energy. In this method, the pressure, p, was 
determined using an explicit relationship of the form p = p(p,e). Consequently, p 
is constant in each element and the shape function x was not used in his case. Also 
a special form of interpolation as given in equation (3.9) was used for internal energy 
in [51]. The finite element system of equations due to Donea et al. [51] is as 
follows : 


Dp = f p (3.11) 

Mv. = f. v (3.12) 
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(3.13) 


E(pe) = f e 


in which each of the entries in the above matrices are given by: 
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In equations (3.10) through (3.13), D and E are diagonal matrices due to the shape 
functions used and the mass matrix M was made to be a diagonal matrix using a 
lumping method [55]. In equations (3.14) through (3.19), denotes an element, 

3 £2 e denotes the boundary of the same element, denotes the applied surface trac- 
tion, and the rest of the notations are the same as before. The dissipation term due 
to molecular viscosity, y(3v^/3Xj + 3v. /3xp, was neglected from the conservation of 

linear momentum equation, equation (3.3). The resulting discrete system of equations 
for this in viscid flow field was integrated in time direction using a predictor- corrector 
method . 

The structural domain was discretized into a number of conical shell elements 
[57] and a time dependent elasto-plasticity finite element method developed by 
Belytschko and Hsieh [57] was used for structural analysis in Donea et al. [51]. The 
structural analysis part will be discussed separately in Section 3.4. 
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One of the example problems considered in Reference 51 was a thin cylinder 
vessel, with hemispherical bottom almost completely filled with water, subjected to 
blast loading from inside. The system of equations were integrated in time up to a 
few milliseconds. The physical problem and the computational results are shown in 
Figure 9 together with moving finite element meshes. Since no viscosity was included 
in Reference 51, a slipping flow boundary condition was used for tangential velocity; 
and the normal velocity was prescribed to be the same as the normal component of 
the moving boundary. This again is obtained as a part of numerical solution at each 
time-level. 

The finite element method for flow analysis used by Belytschko et al. [51,53] 
is basically the same as the one used in Reference 51, but he used an averaging 
method for each of the variables in an element so that the element mass matrix would 
be a diagonal matrix. He also used an hourglass control algorithm to avoid zero- 
energy mode resulting from the averaging method. In this case, the load vector f v 
contains an extra term contributed by the hourglass control algorithm in addition to 
the usual load vectors due to convection, molecular diffusion, surface tensions, and 
body forces. The system of equations was integrated in time- direction , using an 
explicit time integration method, up to a few milliseconds in a finite element modeling 
of hypothetical nuclear reactor core disruptive accident. Their computational results 
compared favorably with the experimental displacement data. In the experiment, 
permanent plastic deformation and severe buckling of steel columns were reported. 

In order to model the buckling of steel columns in the numerical analysis, a slight 
imperfection had been included in the geometry data of the steel column. The 
structural analysis method used can be found in Reference 57. 

Liu solved the same set of flow equations using the penalty method to exclude 
the conservation of mass equation and pressure variable from the system of equations. 
The penalty method has been used very frequently for incompressible laminar flow 
problems [58] and the Stoke's problem [59]. In his method, the stress term, equa- 
tion (3.5), was considered as an extra variable and was interpolated separately, 
which is called the mixed finite element method. He included an elasto- viscoplasticity 
model in the structural analysis part to cover strain rate sensitive solid materials. 

The elasto- viscoplasticity model is covered separately in Section 3.4. 

As can be noticed, there exists significant difference between the nuclear reactor 
safety analysis problems and the flow- solid interaction of a nest of cylinders in cross 
flows. In the latter problem, it is very important to include the turbulent viscosity 
into the flow equations. Any cylinder of a tube bank subjected to cross flows will 
vibrate with respect to its equilibrium position. Even at the onset of failure due to 
large deformation, the distance excursed by any cylinder is small compared to the 
distance traveled by a fluid particle. For these reasons, direct use of the Langrangian- 
Eulerian description may not be necessary for flow- solid interaction analysis of a nest 
of cylinders in cross flows. But the concept of Lagrangian-Eulerian description can 
be used advantageously to obtain convergent solutions for convection dominated flows 
as has been pointed out by Zienkiewicz [58]. 


3.2 Turbulent Flow -Solid Interactions 

Coupled and uncoupled turbulent flow-compliant material interactions were studied 
numerically by Buckingham et al. [16,17] in an effort to reduce drag and aeroacoustic 
noise of submersible hulls. Uncoupled analysis was used to eliminate a number of 
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candidate compliant materials and internal structures (Fig. 10) which do not exhibit 
the desired vibration response from the more elaborate and time consuming fully 
coupled flow-solid interaction analysis. Different materials considered include poly- 
vinylchloride, rubbers, and many others. For all the cases of coupled and uncoupled 
turbulent flow-compliant material interactions considered in References 16 and 17, the 
structural dynamics equation was solved by the finite element codes due to Hallquist 
[60,61]. An elastic or a linear viscoelastic constitutive equation [62] were used to 
model the compliant materials. These constitutive equations and the finite element 
method for structural dynamics are discussed in Section 3.4. 


RIGID SUBSTRATE 



THIN, RELATIVELY STIFF 
NEOPRENE OVERLAYER 

PERIODIC RECTANGULAR 
VOID SLOTS 
SOFT, THICK 
PVC MIDLAYER 


Figure 10. Geometry of compliant surfaces, flow and coordinate system for 
finite element transient response simulation [16,17]. 

In the uncoupled analysis, the fluctuating turbulent wall pressure was generated 
using the Ash code [13,14]. Ash code can generate fluctuating turbulent wall pres- 
sure beneath fully turbulent boundary layer flows. In the method, the instantaneous 
pressure at a spatial and temporal location is obtained by summing up all the pressure 
events affecting the location. Information required to generate each of the pressure 
events include: root mean square value of fluctuating turbulent wall pressure, P rmg 

(or turbulent wall shear stress and its relationship to P rms ) ; probability distribution 

of time between pressure events; probability distribution of frequencies for the 
pressure events; convection speed of the pressure events; decay ratio of amplitude 
of pressure events; wave form of pressure carriers; and etc. In the work of Ash 
[13,14], all this information was obtained from several different sets of experimental 
data for turbulent boundary layer flows; and the statistical results of the generated 
fluctuating turbulent wall pressure compared favorably with experimental data. 
Universality of these experimental correlation equations for different flow situations 
need to be examined further. In order to simulate the fluctuating turbulent wall 
pressure, much statistical information must be known priori for the flow field, there- 
fore his method is of limited use for general flow situations [13]. 

In the fully coupled case [16,17], both the two- and three-dimensional Navier- 
Stokes equations were solved by a pseudo- spectral method [63,64] developed by 
Orszag and Kells [63,64]. The pseudo-spectral method [63,64] was originally 
developed to solve the Navier- Stokes equations for finite amplitude disturbances in 
order to predict transition criterion from laminar to turbulent flow of plane Poiseuille 
and plane Couette flows. The flow field, considered in References 16 and 17, is 
governed by the incompressible Navier-Stokes equations: 


V • v = 0 


(3.20) 
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0 < x < X , 0 < y < Y , and 0 < z i Z (3.22) 

where equation (3.20) is the conservation of mass equation; equation (3.21) is the 
conservation of linear momentum equation; equation (3.22) is the domain of the prob- 

lem;v = ( v ^» v 2 ’ v 3 ^ a column vector of velocities in the three spatial coordinate 

directions; V is the gradient operator; and X, Y, and Z denote extension of the 
domain. The domain in the z- direction was set at a distance of several boundary 
layer thicknesses away from the wall so that the free stream boundary con- 
dition can be used at the outer edge of the domain. The fluid particles at the flow- 
compliant material interface were set to move together in the normal direction, i.e., 
the flow velocity in the normal to the surface of compliant material direction is 
dictated by the normal velocity of the compliant material which is obtained from the 
structural dynamics solution but the flow velocities in both of the tangential to the 
surface directions were set equal to zero. In the free stream flow direction, a 
cyclic boundary condition such that v(0,y,z) = v(X,y,z) was used. In Orszag et al. 
[63,64], the Navier- Stokes equation was first rewritten in a rotational form given as: 
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where u(x,t), u(x,t) = VXv(x,t), is the curl of the velocity vector, and the rest of 
the notations are the same as before. Equation (3.23) was solved by a pseudo- 
spectral method. The pseudo- spectral method has also been used quite frequently 
in the large eddy simulation method which will be discussed in the next section. 
Therefore, a brief explanation and a discussion on the advantages and disadvantages 
of the method is given below. 

The pseudo-spectra method they used belongs to the classical collocation 
method. In the collocation method, solution to the governing differential equations 
are interpolated using basis polynomials and constant coefficients assigned to each 
of the polynomials. In their case, the velocity field was represented by: 


^ ^ " [27r.(mx/X + ny/Y)] 

y = 2^ jL< a ( m * n »P>t) T p (z) e 

|m | <M |n | <N p=o 


where m, n, and p are integers so that a set of (m,n,p) designate each collocation 
point; (2M), (2N), and (P+1) are the number of collocation points in each of the 
spatial coordinate directions respectively; a(m,n,p,t) is function of time only, and 

Tp(z), Tp(z) = cos(pcos ^z) , is the Chebyshev polynomial of degree p. Substituting 

equation (3.24) into equation (3.2.4) and collocating at each of the collocation points 
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yields a system of nonlinear ordinary differential equations, which can be solved by 
any nonlinear ordinary differential equation solution method. But a straight forward 
computation of the collocation method requires enormous computational time as the 
number of collocation points is increased. The special features of the pseudo-spectral 
method due to Orszag et al [63,64] are: orthogonality of the basis polynomials is best 

utilized by performing most of the computations in the Fourier transformed space and 
then mapped back onto the physical space using the fast Fourier Transformation 
method [65] so that the required number of computations are reduced by orders of 
magnitude; and the Adams-Bashforth-Crank-Nicholson time itnegration method devised 
to be used in the pseudo- spectral method is one of the best time integrators for 
ordinary differential equations with regard to truncation errors and numerical stability 
requirement. Detailed explanation on the computational procedure of the pseudo- 
spectral method can be found in a number of literatures [66,67,68]; therefore, it is 
not discussed herein. 

In order to use the pseudo-spectral method to solve any ordinary or partial 
differential equations, basis polynomials which span entire domain and satisfy the 
prescribed boundary conditions must be available a priori. Due to this difficulty, the 
pseudo- spectral method has been used mainly for channel flows or axisymmetric 
annular flows where cyclic boundary conditions can be specified for the upstream and 
the downstream locations. The validity of cyclic boundary conditions has been 
questioned in References 67, 68, and 80. The boundary conditions are believed to 
be approximately applicable for fully developed channel or pipe flows [63,64,92]. 
Accordingly, the method may be suitable for problems with extremely limiting cases of 
boundary conditions defined on simple regular domains. Application of the method 
to flow-solid interaction analysis of a nest of cylinders subjected to cross flows would 
be extremely difficult, if not impossible, due to the complicated domains and boundary 
conditions . 


3 . 3 Turbulence Models 

In numerical analysis of a nest of cylinders subjected to cross flows, 
pressure straining and redistribution of Reynolds stresses near the wall region may 
become important due to a number of cylinders in the cross flow. Also in flow-solid 
interactions, fluid force is transmitted to the structure through wall pressure and 
skin friction; therefore, evaluation of fluctuating turbulent wall pressure is very 
important. Since successful numerical analysis of flow-solid interaction of tube banks 
would require an appropriate unsteady turbulence model to be incorporated into the 
flow analysis, a few turbulence models are discussed below. 

Significant progress has been made in modeling steady turbulent flows. But due 
to pressure straining, redistribution of Reynolds stresses, and the complicated geome- 
try of tube banks, some of the steady state turbulence models with proven performance 
may not even be applicable for steady turbulent flows through a nest of rigid cylin- 
ders. For example, low Reynolds number turbulence models which include the dis- 
tance from the walls explicitly in the turbulence equations may hardly be applicable 
for flows through tube banks, since there exist so many cylinder walls in any tube 
bank, e.g., there are 600 tubes in the SSME Main Injector Assembly. On the other 
hand, numerical modeling of unsteady turbulent flows has shorter history than the 
steady case; consequently, there exist only a limited number of publications in this 
area and physical understanding of unsteady turbulent flows is not complete at the 
present time. Therefore, this discussion is limited to the most successful turbulence 
models which may be further extended to unsteady turbulent flows and which can be 
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used for turbulent flow-solid interaction of tube banks. These are a few low Reynolds 
number k- e turbulence models, a few different versions of the algebraic Reynolds 
stress models, and the large eddy simulation method (LES). The LES method cannot 
be used for turbulent flows through tube banks at its present state of development as 
can be seen in the following discussions, but the method has been originally developed 
for and primarily used for unsteady turbulent flows. Comprehensive reviews of 
various turbulence models can be found in Nallasamy [69], Martin [70], Lakshminarayana 
[71], and the references appearing in these literatures. 


3.3.1. Extended k- e Turbulence Models and Algebraic Reynolds Stress Models 

For a class of complex flow problems, anisotropy of Reynolds stresses induce 
significant differences in the predicted flow field, e.g., the existence of secondary 
flow in curved square ducts [72]. For these flows, the conventional (or standard) 
k- e turbulence model yields unsatisfactory computational results. The Reynolds stress 
model [73] can describe the convection, diffusion, and anisotropy of each component 
of the Reynolds stresses. But use of the Reynolds stress model is prohibited in many 
cases of practical importance due to the number of variables which result from six 
differential equations for Reynolds stresses and one equation for turbulence dissipation 
function for three-dimensional flows. 

In order to account for the anisotrophy of Reynolds stresses, the algebraic 
stress model, which retains most of the basic features of the original Reynolds stress 
model, was developed [74]. The algebraic relationship between the six Reynolds stress 
components are obtained from the Reynolds stress equation by neglecting the convec- 
tion and diffusion terms [76] or by assuming the convective transport of Reynolds 
stresses to be proportional to the transport rate of turbulent kinetic energy [74]. 

Hence in algebraic Reynolds stress models, the k-e turbulence models are enhanced 
with this algebraic relationship instead of using the classical Boussinesq eddy viscosity 
assumption [76]. 

In the numberical analysis of turbulent flows, the wall boundary conditions are 
usually applied in two different ways. In the first approach, the computational 
boundary of flow domain near the wall is slightly separated from the wall and the 
wall functions are applied at the first grid point. Strictly speaking, these wall 
functions are not applicable near the separation points on the walls; nevertheless, 
they have been used successfully in many cases. In the second approach, the first 
grid point is located on the wall, no slip boundary condition is applied at the grid 
point, and the necessary modifications are included in the turbulence models in order 
to account for the low Reynolds number effect near the wall. 

In flow-solid interaction of a nest of cylinders subjected to cross flows, the 
vibration of cylinders is a three-dimensional motion. The capability of wall function 
methods for the analysis of flows with moving boundaries has not been established as 
yet, and need to be established so that the method can be used for the analysis of 
flow-solid interactions. Alternatively, the low Reynolds number turbulence model may 
be used so that the no-slip flow boundary conditions can be specified at the flow- 
solid interfaces. 

For unsteady turbulent flows, the mean and fluctuating quantities are defined 
differently than those for steady turbulent flows. Therefore, it would be appropriate 
to include this topic. The most familiar averaging technique for unsteady turbulent 
flows would be the ensemble averaging method in the sense that it yields exactly the 
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same flow equations as the time averaging method for steady turbulent flows. Ensemble 
average of any fluctuating quantity, V(x,t), is defined as: 

N 

¥(x,t) =^ 1 ™ oo ^ f(x + nx) (3.25) 

n=0 


where x is the period of fluctuation, n is an integer, and V denotes ensemble average 
of T . On the other hand , the time average for steady turbulent flows is defined as 


t+T 

?(x,t)=^ oo i f Y(x,t) dt 

t 


(3.26) 


where V denotes the averaged quantity. For both of the averaging techniques, the 
velocity component is decomposed as 


y = X + y’ (3.27) 

T 

where v = ( v 2 » v 2 ,v 3 ^ * s a co ^ umn vector of the three velocity components in the 

three spatial coordinate directions, and v’ is the corresponding column vector of the 
fluctuating velocities. The time dependent turbulent flow equations for incompressible 
flows are given as: 
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= 0 (3.28) 
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(3.29) 


The Reynolds stress, v.’v.’, in equation (3.29) is modeled using the Boussinesq 
viscosity assumption given as: 


•Vi'Vj’ 
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6 i) k 


(3.30) 


where is the turbulent viscosity; 6- is the Kronecker delta such that <5^ = 1 for 
i = j, and = 0 for i ^ j; and k is the turbulent kinetic energy, k = 1/2 V.'v.' . 
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The standard two equation (k-e) turbulence model due to Launder and Spalding 
[78] is given as: 


ak — ak _ a 

at v j ax. ax. 
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(3.31) 


(3.32) 


where e is the dissipation rate of the turbulent kinetic energy , 


e = v 


av.* av.' 


ax. ax. 


P is the production rate of the turbulent kinetic energy, 


P = 



av. 
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ax. 
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V. = c 
t t VI 

= 0.09, c 


k /e, is the same turbulent viscosity as used in equation (3.20); and 
= 1.44, c 0 = 1.92, = 1.0, 


. . - x.t. 1 , v. „ = 1.92, a, = 1.0, and a =1.3 have been used in Reference 

Vi 'el e2 k e 

78. Since equations (3.31) and (3.32) are not valid in the viscous sub-layer region 
near the wall, where the molecular viscosity is dominant over the eddy viscosity, the 
wall function methods are used to supply wall boundary conditions for the differential 
equations [78]. There exist a few modified k-e turbulence models to include low 
Reynolds number effects. The extended k-e turbulence model due to Hassid and 
Poreh [79] is given as 
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(3.34) 


28 


where v, = c k 2 /e, c =c f , f =1- exp (-0. 0015R. ) , R. = k 2 /(ve) is the tur- 
t y ’ y y °° y y ^ v t' ’ t 

bulent Reynolds number, and c ^ = 0.09. The model by Jones and Launder [80] is 

given as y 



(3.36) 


where ^expf-2. 5/( l+R t /50)] , c ^ = C £oo [l-O.3exp(-R t 2 )] , c = 0.09, c gl = 

1.45, and c ^ = 2.0, and rest of the notations are the same as in equations (3.33) 

and (3.34). The wall boundary conditions for the low Reynolds number turbulence 
equations are given as: 


k = 0 
e = 0 


(3.37) 


Various algebraic relationships between Reynolds stresses can be obtained by 
introducing simplifying assumptions into the Reynolds stress equation [73,77]. These 
simplifying assumptions are: production and dissipation of turbulent kinetic energy 

are in equilibrium, transport of Reynolds stresses are proportional to the transport 
of turbulent kinetic energy, transport and diffusion of Reynolds stress are negligible, 
and any appropriate combination of these assumptions. The Reynolds stress equation 
[73,77] and the algebraic Reynolds stress model due to Rodi [74] are given below. 

The Reynolds stress equation [73,77] is given as: 
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(3.38) 


where (I) is the time rate of change, (II) is the convective transport, (III) is the 
stress production, (IV) is the pressure straining term, (V) is the diffusive transport, 
and (VI) is the dissipation of the Reynolds stresses. The closure model for higher 
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order terms for low Reynolds number flows proposed by Hanjalic and Launder [77] 
are given as 
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where P = -v.'v, ' (9v /9x. ) is the rate of generation of turbulence energy by the 

36 K 36 1 o 

mean strains; f = (1+R./10) ; R = k7(ve) is the turbulent Reynolds number; and 

S X X 

c^ = 1.5, C 2 = 0.4, and c g = 0.11 have been used in Reference 77. The turbulence 
dissipation function, e, in equations (3.39) through (3.45) is determined by equation 
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(3.46) which is the convection -diffusion equation of the turbulence dissipation func- 
tion. The pressure straining due to wall proximity, equation (3.40), was proposed in 
Reference 73 but was disregarded in Reference 77. The form of wall proximity equa- 
tion, equation (3.40), is not appropriate for numerical analysis of turbulent flows 
through the tube banks due to the explicit appearance of the distance from the wall. 

The algebraic Reynolds stress model derived by Rodi [74] is given below, where 
the convective transport term of the Reynolds stresses is assumed to be proportional 
to the transport of turbulent kinetic energy, and the diffusion of Reynolds stresses 
is neglected. Accordingly, 


Vi’Vj’ 
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(3.47) 


where y % 0.6, and the rest of the notations are the same as in the Reynolds stress 
equation. Successful application of the algebraic Reynolds stress models can be found 
in a number of publications [74,81,82]; hence, it is not elaborated herein. But its 
extension to unsteady turbulent flows remains to be tested. Computational time for 
using algebraic Reynolds stress models is comparable to the k-e turbulence models 
for most of the cases. 

Extension of the k-e turbulence models to unsteady flows has been limited to 
flows with simple geometries [83] and boundary layer flows. Even for these cases, 
straight forward extension of the k-e turbulence models has been more unsuccessful 
than successful. These model inadequacies motivated a group of researchers to devise 
a new triple decomposition model of turbulence [ 84] , in which any fluctuating quantity 
is decomposed into three parts, i.e., time mean quantity, organized oscillations, and 
random fluctuations. However, this method also remains to be further studied to be 
useful for engineering calculations. 

A successful application of the standard k- e type turbulence model to unsteady 
recirculating turbulent flow can be found in Boyle and Golay [85]. Measured velocity 
field and computed velocity field at 80 sec into transient of the scaled nuclear reactor 
outlet plenum are shown in Figure 11 [87]. Experimental and computational 
velocity components, at various locations inside of the plenum chamber, which vary in 
time are shown in Figure 12 [87]. The computational results compare favorably with 
experimental data during the flow development period, but the same results begin to 
diverge from experimental data as steady state is approached. 


3.3.2 Large Eddy Simulation (LES) Method 

In turbulent flows, instability of mean flow generates large eddies and these 
large eddies break up into smaller eddies. After many cascade processes, see Figure 
13 [86], the fine scale eddies are dissipated by viscous forces. The scale of large 
scale turbulence, or motion of large eddies, strongly depends on the type of flow and 
geometry of the flow field; whereas, fine scale turbulence, or motion of fine scale 
eddies, is almost universal for any turbulent flows. Therefore, fine scale turbulence 
can be modeled more easily than the large scale turbulence due to its universality. 

The large eddy simulation (LES) method best utilizes these physical observations. 

Until now, the LES method has been used for channel flows and axisymmetric annular 
flows. 
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Figure 11. Turbulent recirculating flow in a nuclear reactor outlet plenum [85]. 
(a) Measure mean velocity profile after 80 sec into transient, (b) Computed 
mean velocity profile after 80 sec into transient. 
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Figure 12. Comparison of measured vertical velocity with computer 
velocity [85], a through h in the Figure denotes spatial 
locations specified in Figure 10(a) . 
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Figure 13. The energy cascade process [86], 
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In the large eddy simulation method, any flow variable is averaged over a finite 
volume which is determined by the grid size used in computations. This process is 
called filtering. Let ¥ be some flow variable, then Y can be decomposed as: 


Y = Y + Y' 


(3.48) 


where Y’ is called the residual field, and Y is the large scale field defined as: 


nx,t) 
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AV 


G ( x , x ' ) Y (x' ,t) dx 


(3.49) 


where AV = Ax^AXgAXg, Ax k (k=l,2,3) is the grid size in the k-th spatial coordinate 

direction, and the integration is performed over the volume AV. The filtering func- 
tion, G(x,x'), in equation (3.49), was first introduced by Leonard [87] in order to 
account for nonnegligible variation of large scale motion inside of the averaging volume, 
AV. The filtering function is given as: 


G(x,x') 


G^(x^,x^') Gg(x2 ,X2') G^x^jX^ ^ 


(3.50) 


where Gj(i=l,2,3) is the distribution function. The most frequently used distribution 
function is the Gaussian distribution function given as: 
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(3.51) 


where A is twice the computational mesh size in x^- coordinate direction. Many differ- 
ent forms of filtering functions can be used in practice. For example, G(jc.x') = 1 
was used in the first application of the LES method in Deardroff [88]. A flux form of 
the Navier- Stokes equation used in LES is given as [88]: 
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Filtering equation (3.52), using any filtering function G(x,x'), yields: 
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(3.53) 
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Kronecker delta. A s ubgrid- scale-model (closure model for LES equation) of the 

residual stress, v.’v .' , due to Deardroff [88] is given by: 
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where v.v.’ = v.v.' = 0, and c. is a constant to be determined experimentally. Lilly 


i 1 


1 i 


[89] proposed c^ = 0.17, Deardorff [88] used c = 0.1, and Moin and Kim [90] used 

c^ = 0.065. In equation (3.53), n 6^/3 has been added to the pressure gradient 

term to compensate for the same term appearing in the expression so that equation 

(3.55) holds when the subscripts i and j are contracted. Consequently ru . in equa- 
tion (3.53) need to be modeled separately, and Lilly [89] proposed: 
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where c^ = 0.094. 


The governing equation for pressure distribution 


By letting P = P/P + n kk 
tion (3.53), yields: 
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can be obtained as follows, 
the momentum equation , equa- 
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where dissipation due to molecular viscosity has been neglected and the time derivative 
term can be dropped out using the conservation of mass condition. 
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Several different methods have been used by different workers [88,91,92] to 
solve the LES equation, and the form of the LES equations were slightly modified to 
be suitable for each of the numerical methods used. In Deardorff [88], the classical 
finite difference method was used to solve the LES equation as given in equation 
(3.53); in Schumann [91], the control volume finite difference method [90] was used 
and any time derivative term appearing in the Navier- Stokes equations was filtered 
over the control volume as before whereas the rest of the terms appearing as fluxes 
in the control volume finite difference method was filtered over each of the surfaces 
of the control volume where fluxes are evaluated. In Moin and Kim [92], the pseudo- 
spectral method [63,64] and the LES equation of the rotational form of the Navier- 
Stokes equations, i.e. , equation (3.23), were used. Particularly in Moin and Kim [92], 
the LES equation was refined to include the exact treatment of the Leonard’s stress, 
i.e., equation (3.54), and molecular viscosity. 

The computational results given in these references [86,88,92] compare 
excellently with experimental data. But it should be noted that the constant coeffi- 
cient c^ , used in equation (3.56), as proposed by Lilly [89], is larger than twice the 

value used in Moin and Kim [92]. Hence, the eddy diffusity, which is proportional 
to the square of c A , used in these calculations for the similar channel flows is quite 

different. Conceptually, the subgrid-scale-model should not dominate the computa- 
tional results when fine grids are used for numerical analysis, and the residual stress 
itself should vanish as computational grid becomes infinitely small. But the afore- 
mentioned computational results, using practical grid size, showed that the subgrid- 
scale model is important to render the computational results comparable to experi- 
mental data. Therefore, further refinement of the sub grid -scale -model need to be 
made for the LES method to be applicable to practical engineering calculations. 


3.4 Structural Analysis Methods Including Large and/or Permant Deformations 

Due to fluidelastic instability, thin cylinders in tube banks may go through 
large deformations which may cause viscoplastic deformations of the cylinders. There- 
fore, the structural analysis method to be used for flow-solid interaction of tube banks 
need to be capable of modeling elastic vibrations as well as large and/or permanent 
viscoplastic deformations, in order to predict Bauschinger effect, hysteretic effect, 
and cyclic straining. But to date, no specific numerical structural analysis method 
designed for this purpose has been developed. Structural analysis methods which 
can be further developed to meet these requirements are discussed in this section. 

Most of the structural analysis methods are based on the finite element displace- 
ment method. The finite element system of equations are derived through the mini- 
mization of energy functional or derived from the virtual work principle. Elasticity 
problems are usually described on the frame work of Eulerian coordinate system; and 
hyperelasticity, large deformations, and viscoplastic deformation problems are usually 
described on the convected coordinate system in order to account for large displace- 
ments. In the methods utilizing the convected coordinate system, the virtual work 
principle is described on the Lagrangian coordinate system which is updated at each 
of the solution time levels. The virtual work principle and the minimization of energy 
functionals are identical for elasticity problems; but the virtual work principle is more 
convenient for geometrically and/or materially nonlinear problems, since derivation of 
energy functionals for nonlinear problems can be avoided. In the following discussions, 
the virtual work principle and the convected coordinate system are used. 
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For time dependent deformation of structures, consider a point initially at 
x.(i=l,2,3) moved to u.(l=l,2,3) in the same cartesian coordinate system. In Lagran- 

gian formulation, the deformations are described in terms of the convected coordinates 
Uj(i=l,2,3). The force equilibrium equation, or the momentum equation, for a struc- 
tural system is given as 


a.. . + b. = 0 
1M i 


u(x,t) = u Q (t) on 


cr.. n. = T.(t) on 3£2 9 
ij ] i * 


(3.59) 

(3.60) 

(3.61) 


where bj is the body force; Q is the structural domain; 3^ is the boundary on which 
the Dirichlet boundary condition is specified; 3 fig is the boundary on which surface 
traction is specified; 3fi^f") 3^ * s a nu ^ se t> 3 fig = 9ft, and 3 fi is the entire 

boundary of the domain; both the surface traction T.(t) and the prescribed displace- 
ments Uq( 0 are functions of time in general; and n is an outward normal vector on 

the boundary. The inertial force due to acceleration is usually included into the 
momentum equation as a negative body force. The virtual work statement is given as: 
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where 6u.(i=l,3) is the virtual displacement which satisfies the Dirichlet boundary 

condition, equation (3.60), on the boundary 3^. Integrating by parts the stress 
term and using the fact that 
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a.jn.iSUjdS = 0 , 


since 6u. = 0 on 3fi^ in equation (3.61), yields: 
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(3.62) 


which is the virtual work expression to be used in the finite element method. In 
finite element analysis, it is more convenient to write the stress tensor, a-, and the 
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strain tensor, e.., into column vectors c T = < a xx > CT yy > CT ZZ > o xy >° yz > a zx >» and e, T = 

{e xx> e yy ,e zz ,Y xy ,Y yz ,Y zx }. The shearing strain tensor, e..(#j), is not equal to the 

engineering shearing strain, hence adequate modification need to be included 

in the material property coefficients in the stress-strain relationship. Using the 

T 

matrix notations, o..6Uj ^ can be expressed as 6e g. Let 
u = N a 

i — py r-* 

e = L u (3.63) 

<— 

e = B a 

where N is a matrix of interpolating polynomials, a is the column vector of nodal dis- 
placements, L is a differential operator, and B = . For geometrically nonlinear 

problems, L and B depend on displacements. Introducing equations (3.63) into equa- 
tion (3.62) and taking the variation of equation (3.62) with respect to the nodal dis- 
placements yields , in matrix notation , 



dx - 
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N 1 b dx - 



T dS = 0 


(3.64) 


For geometrically and/or materially nonlinear problems, equation (3.64) is not 
generally satisfied at any computational stage. Hence, writing equation (3.64) in 
incremental form yields 
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B a dx + 


f B T da + f N' 
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(3.65) 


where dT = 0 if the surface traction is not incremented at any computational stage. 

Equation (3.62), or equation (3.65), represents the most general form of the 
finite element equation for any spatial dimensions; however, its specific forms for 
different problems and different solution methods are so diverse that they cannot be 
covered in this review. These different problems include: different strain/stress 

assumptions; geometrical differences such as one dimensional bars, two dimensional 
plates, and quasi three dimensional shells, etc.; and the geometrical and/or material 
nonlinearities. Different solution methods include: interpolation methods such as 

Lagrangian interpolating polynomials and the Hermite interpolating polynomials; 
different formulations such as mixed formulation and displacement formulation; and 
even time integration methods. Only a couple of algorithms that can be extended 
easily for the flow-solid interaction analysis are covered herein. These are the 
elasto- viscoplasticity analysis method of Hughes and Taylor [93] and the nonlinear 
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vibration analysis method, which includes large deformation and elasto- viscoplastic 
deformations, due to Rodal et al. [99,100]. 

The most rigorous mathematical treatment of viscoplastic deformation can be 
found in Hughes and Taylor [93] and Owen and Hinton [94], In References 93 and 
94, a linear strain- displacement relationship has been used and the acceleration force 
has been neglected, on the assumption that it is negligible compared with forces which 
cause plastic deformation. Hence, the method is suitable for problems in which only 
the permanent plastic deformation is important, but not suitable for problems in which 
the vibrational aspect is as important as the permanent plastic deformations. The 
general form numerical analysis method for viscoplasticity due to Hughes and Taylor 
[93] and Owen and Hinton [94] is given below. 

In elasto- viscoplastic analysis, it is assumed that the total strain rate, e, can 
be decomposed into an elastic part, e^, and a viscoplastic part, £ V p> i.e.. 


£ = 


£ + £ 

rJ vp 


(3.66) 


where the over-dot, (‘) notation denotes time rate of change. The viscoplastic strain 
rate due to Perzyna [95] is given as: 


e = y<cf>(F)> 

—vp ' Y v ' 3 a 


(3.67) 


where y is the fluidity of the plastic flow; <<}>(F)> implies < 4> (F) > = <J> (F) if F > 0, and 
<<HF)>=0ifF <0;Fis the yield function such as the Huber- von Mises yield func- 
tion [55] , and Q is the plastic potential function. For associated plasticity, F is 
identical to Q , but different for nonassociated plasticity . The stress rate is related 
to the elastic strain rate as: 


a = D i = p (£ - e ) (3.68) 

or, in incremental form, 

Aa n = D(Ae n - Ae" ) (3.69) 
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where the superscript n denotes the time level t = t , A denotes increment, and the 

increment is defined for a time interval from t tot , i.e., At = t ,,-t . The 
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viscoplastic strain increment , A e 11 , can be interpolated as 
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where 6 is an interpolation parameter such that 0 = 0 corresponds to the fully explicit 
time integration method, 0 = 1 corresponds to the fully implicit method, and 6=1/2 
corresponds to the Crank -Nicholson type method. The time integration method for the 
present elasto- viscoplasticity analysis is unconditionally stable for 0 > 1/2. As a 
remark, numerical stability does not imply that convergent solutions can always be 
obtained for many nonlinear problems. Since the viscoplastic strain rate is a function 
of stress state, equation (3.70) can be rewritten as 


. n . , 
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where (e”p 1- e” p ) = C 3 £ V p / 3 A 2^ ^ as been used, and c n = 0At n O e vp /9 a) n . The 
incremental strain-incremental displacement relationship can be written as: 


Ae n = B Aa 11 (3.72) 

/V- 4 ' 

where B is a constant matrix for infinitesimal strain- displacement relationship. Sub- 
stituting equations (3.71) and (3.72) into equation (3.69) yields: 


Aa n = D n (B Aa n - *e n At ) (3.73) 

^ ri n — 1 

where D = (D + C ) Equation (3.73) is the incremental constitutive equation 
to be used in~the incremental virtual work principle, i.e., equation (3.65). For 
material nonlinearity only , the incremental virtual work equation is given as : 
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Inserting equation (3.73) into equation (3.74) yields 
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The incremental displacement, A a is obtained by matrix inversion of equation (3.75) 
as 


A a 


n 


where 

^ X 



= B T D n B is called the 


At dx - 
n r- 

tangential 


Af" 

stiffness matrix. 


(3.76) 
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In Hughes and Taylor [93], the above methodology was used to solve an elasto- 
viscoplastic deformation of a thick wall cylinder subjected to constant internal pressure. 
The physical problem was reduced to a one- dimensional problem by assuming an 
axisymmetric plane strain state. The physical problem and the computational results 
are shown in Figure 14. 



THICK-WALLED CYLINDER 


^ 4—4 

0.125 

(a) 



a = 1.0 1 -= t = 10 6 


RESULTS FOR TIME STEPS 
2-12 ARE IDENTICAL 


(b) 


Figure 14. An elasto-viscoplasticity analysis of a thick cylinder subjected 
to internal pressure [93]. (a) Description of physical problem and 

a finite element discretization, (b) Comparison with 
exact plasticity solution. 
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The numerical analysis method given in equations (3.66) through (3.76) are so 
general that it can be used for any elasto-viscoplastie problems. Yet, specific forms 
of strain-displacement relationships for individual problems need to be derived. For 
flow- solid interaction of tube banks, a general shell element would be preferred. 

Refinement of the lasto-viscoplasticity analysis methods in References [93] and 
[94] can be made in several different ways. These refinements would be: incorpora- 

tion of temperature dependency of the viscoplastic strain rate into equation (3.67); 
inclusion of geometrical nonlinearity; inclusion of work hardening into the yield func- 
tion; and inclusion of the acceleration force term into the governing equation so that 
Bauschinger effect, hysteretic effect, and cyclic strain hardening can be predicted. 
Examples of temperature dependent viscoplastic strain rate can be found in References 
96 and 97, in which flow formulation was used and the effect of elasticity was com- 
pletely neglected. Examples of geometrical nonlinearity and work hardening can be 
found in References 98, 99, and elsewhere. 

Rodal et al. [99,100] developed a finite element structural analysis method which 
includes finite strain, large deformation, and elasto-viscoplasticity so that Bauschinger 
effect, cyclic straining, and hysteretic effect could be predicted. In this method, all 
the nonlinearities such as geometrical and material nonlinearities were incorporated into 
the load vector term so that the contributions of the nonlinearities appear only as 
modifying forces to the applied load. The structural dynamics equation due to Rodal 
et al. [99,100] is given as 

M a + K a = f + f 11 *' (3.77) 


where M is the usual mass matrix, K is the elastic stiffness matrix, f is a column 

vector of externally applied loads, and f 11 ^ is a column vector of modifying loads due 
to geometrical and material nonlinearities. The computational procedures used in 
References 99 and 100 are simpler than the one used in References 93 or 94, in which 
the viscoplasticity has been treated more rigorously . Undoubtedly , the numerical 
method due to Rodal et al. [99,100] would be adequate for modeling short period 
vibrations such as those associated with impact problems; but the accuracy of the 
method as applied to structures vibrating for a longer period of time than the impact 
duration need to be tested further. 

The time dependent elasto-plastic finite element method [57] used in the laminar 
flow- solid interaction analysis of Donea et al. [51] and Belytschko et al [52,53] are 
described below together with the Newmark time integration method [101]. The 
implicit Newmark time integration method is suitable to achieve a rigorous force equi- 
librium condition in integrating a system of second order ordinary differential equa- 
tions in the time coordinate direction. But its adequacy for flow-solid interaction has 
to be tested yet. Elasto-plasticity constitutive equations can be found in a number of 
publications including References 55 and 102, and are not included herein. 

Belytschko and Hsieh [57] developed a time dependent elasto-plasticity finite 
element method. In this method, infinitesimal strain assumption (a linear strain- 
displacement relationship) was used and the large dispalcement was taken care of by 
using the convected coordinate system. The finite element equation used in Reference 
57 is given as: 
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(3.78) 
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p w u dx + 
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a dx 



T dS 


which is equivalent to equation (3.64), but the acceleration term appears explicitly. 
The semi- discretized (fully discretized in spatial coordinates, but left continuous with 
respect to the time variable) finite element equation becomes 


M a 


- K 



(3.79) 


In Belytschko and Hsieh [57], the mass matrix was lumped into a diagonal matrix and 
the elasto-plastic stiffness term was incorporated into the load vector term as a modi- 
fying load vector. The resulting system of second order ordinary differential equa- 
tions was solved explicitly in time. As discussed before, explicit time integration 
methods are suitable to solve short-period structural responses such as the one for 
impact problems, or to solve transient problems which asymptotically approach steady 
state. The force equilibrium condition may not be satisfied rigorously at each time 
step. The implicit Newmark time integration method which can satisfy force equilibrium 
condition more rigorously is briefly described below. 

[101], the displacement vector in equa- 
time, is expanded into a parametric 


(3.80) 


where 

c Q = l/[y(At n ) 2 ] 
Cj. = l/(y At n ) 
c 2 = l/2y - 1 


In the Newmark time integration method 
tion (3.79), which is a continuous function of 
series form such that 


•n+1 _ , n+1 n, n -n 

^ c Q (a ^ ) c x a c 2 a 


a n+1 = a" + [(1 - 5) 5 n + Sa n+1 I at 


n 


a n+1 = a n + a n At + [(0.5 - y) ‘a 11 + y 


and a and y are constant such that (a,y) = (1/2, 1/4) corresponds to the uncondi- 
tionally stable, constant- average-acceleration time integration method. Substituting 
equations (3.80) into equation (3.79) yields 
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(3.81) 


(c 0 M + g) a n+1 = c Q M [a n + (c^Cq) a n + (c 2 /c Q ) a 11 ] + f 


which can be solved for displacements for each new time-level. 


3.5 Numerical Solutions of Flows Through a Nest of Cylinders 

Only very recently , a group of researchers began to numerically model the fluid 
flow through tube banks. The recent developments of flow analysis method for tube 
banks are discussed in this section. These are: control volume finite difference 

analyses due to Singhal et al. [44], Wambsganss et al. [103], and Sha and Yang [104]; 
two-dimensional finite element analysis due to Zienkiewicz et al. [105]; and body fitted 
grid finite difference analyses due to Kwak and Chang [106], Chang and Kwak 
[107], and Rogers et al. [108]. In all of the above studies, the flow fields were 
assumed to be steady and the viscosity is constant. 


3.5.1 Control Volume Finite Difference Analyses of Steady Laminar Flows Through 
Tube Banks 

For tube banks with a large number of cylinders, detailed numerical modeling of 
the flow field is prohibited due to the present computer capability including available 
memory and computational speed. Approximate numerical analysis methods based on 
the concept of flow through porous media were proposed by Sha et al. [104], 
Wambsganss et al. [103], and Singhal et al. [44]. In these methods, the flow domain 
is discretized into a number of cells where each cell contains several cylinders (Fig. 
15) . Then the resistance of each cell to fluid flow is represented by volume porosity 
(fraction of volume occupied by the fluid in the cell) and surface permeability (frac- 
tion of open projected flow area in the direction of fluid flow in the cell) . The 
governing equations are the same Navier- Stokes equations, but any volumetric quan- 
quantity needs to be multiplied by the volume porosity; and each velocity component, 
by the surface permeability in the direction of each velocity component. The flow 
field shown in Figure 8 [44] was obtained by solving the flow equations for porous 
media using the control volume finite difference method [90]. 

In the fluidelastic instability analysis of Wambsganss et al. [103], the computed 
flow velocity normal to the cylinder walls was taken as reference velocity for predict- 
ing the fluidelastic instability using the instability diagram shown in Figure 4. The 
validity of these methods remains in question; since fluidelastic instability is caused 
primarily by flow-solid interaction, and in most of the cases, the most serious fluid- 
elastic instability is confined to the flow inlet region of tube banks according to 
experimental observations. 

More recently, a new computational method, called the control- volume based 
finite elements, was developed by Baliza and Patankar [116], Prakash and Patankar 
[117], and Prakash [118] among many others. The method utilizes isoparametric 
finite elements to discretize the domain and to interpolate the solution vectors. The 
discrete system of equations is derived by considering the flux through the element 
sides as in the control- volume based finite difference method. On the other hand, in 
the standard finite element methods, the discrete system of equations are derived 
using the method of weighted residuals (MWR) to be introduced in the following sec- 
tion. The mathematical foundations of the two different methods are quite different; 
and consequently, the final discrete systems of equations are also quite different. 
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Figure 15. Control volume finite difference analysis of flow through 
a porous media [103]. (a) Discretization of a tube bank. 

(b) Configuration of a cell. 

Until now, the control-volume based finite element method utilized only the linear iso- 
parametric elements. In this case, the curved geometry is approximated by straight 
lines connecting the nodal points. On the other hand, most of the finite element 
methods for structural analysis utilizes quadratic or higher order elements to best 
approximate the geometry as well as solution vectors. Use of the control -volume based 
finite element method for flow -solid interaction analysis in the future remains as an 
open question. 
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3.5.2 Finite Element Analysis of Steady Two-Dimensional Laminar Flows Through 
Tube Banks 

The finite element method has certain advantages over the finite difference 
method. These advantages are: capability to include various types of boundary con- 

ditions more naturally and rigorously in the formulation; capability to show convergent 
nature of the method mathematically; and capability to model complicated domains more 
precisely than the finite difference method. 

Application of the finite element method for flows has shorter history than the 
finite difference method, and the finite element method is still in its early development 
stage as far as application to fluid flows is concerned. Much of the development time 
for the finite element method for flows has been spent in overcoming numerical insta- 
bility arising from convection in high Reynolds number flows. But with the recent 
development of the Petrov- Galerkin method [109] for convection dominated flows, more 
rapid development of the method for flows is expected in the near future. The 
Petrov- Galerkin method and its application to two-dimensional steady laminar flows 
through tube banks are discussed in this section. 

The Bubnov -Galerkin method [109] (in which test functions are selected from 
the same space as the trial functions) yields exactly the same discrete equation as 
the central differencing finite difference method for one- dimensional problems discre- 
tized using linear elements. As a remark, for most of the cases, the discrete equa- 
tions obtained through the two different methods are quite different. For a long time, 
it was believed that the Bubnov- Galerkin method (or equivalently, the central differ- 
encing scheme) is the most accurate discrete representation of the continuous differ- 
ential equations. But, Brooks and Hughes [109] showed that: the Bubnov- Galerkin 

method underestimates the numerical diffusion; the conventional Petrov-Galerkin method 
(or equivalently, the backward differencing scheme) overestimates the numerical 
diffusion; and the optimal Petrov-Galerkin Method (this should be considered as the 
consistent Petrov-Galerkin method when source terms exist in the differential equation) 
yields the most accurate numerical results. A finite element solution of one- dimensional 
convection -diffusion equation of the form: 


r. rp J\ 

V yr " k 0 TT = ° in Xe(0 ’ 1) 

3x 


(3.82) 


with boundary conditions T(0) = 0 and T(l) = 1, is shown in Figure 16, where T is 
the temperature. The basis of test functions for the consistent Petrov-Galerkin method 
is given as [105] 
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(3.83) 


kg" = a v h/m 


(3.84) 
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Figure 16. A finite element solution of one -dimensional 
convection -diffusion problem. 

where a = coth(P e /2)-2/P e , P g = |v|h/k is the grid Peclet number, v is the convec- 
tion velocity, <Jk(x) is the shape function, m = 2 for linear elements, and is the 

number of nodes in an element. This method naturally extends to two-dimensional 
flows as shown in the two-dimensional flows through tube banks. Extension of the 
method to three-dimensional cases was suggested to be the same [105]. 

For a two-dimensional case, the basis of the test fucntions of the consistent 
Petrov-Galerkin method is given as [ 105] : 


^(x) 


<f>i*(x) 


3+i 



l = 1,2 


(3.85) 


T 221/2 

where x = (x,y), k^" is the same as in equation (3.84) but |v| = (v^ +Vg ) at 

each of the Gauss points used in numerical integration, and h is the length scale of 
the element . 

The finite element equation for Navier- Stokes equations is obtained by the 
method of weighted residuals. In the method: the flow domain is discretized into a 

number of elements; the Navier-Stokes equations are multiplied by the test functiQn; 
the second order diffusion term is integrated by parts; continuous variables are inter- 
polated using the nodal values of the variables and interpolating polynomials; the 
resulting equations are integrated over each element using the Gauss numerical quad- 
rature to obtain element system of equations ; and the element system of equations are 
assembled together to obtain a global system of equations. A standard finite element 
system of equations for an element is given in the following. 
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integration is performed over each element, and <)>. x denotes 3<}>j/3x, etc. 

In practice, the penalty finite element method [119] has been used quite fre- 
quently for flow problems (or more frequently for some problems such as polymer 
extrusion problems [110]). In the penalty method, the unknown pressure is pre- 
eliminated from the Navier- Stokes equations. The pressure field solution can be 
obtained in the post process if necessary. The laminar flow through tube banks 
shown in this section was also solved using the penalty method. To eliminate the 
pressure variable , we set : 


V • v = P/< 


(3.89) 
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where p is the pressure, and <*p is the penalty parameter. As ap approaches infinity, 

X • v ^ 0 and, hence, the continuity condition is satisfied approximately. In these 

7 10 

computations, is taken to be 10 through 10 times the molecular viscosity to avoid 

ill-conditioned system of equations. Finite element discretization of equation (3.89) 
yields 


[Q 


T 
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(3.90) 


where and Q2 are the same as in equation (3.88), S-J Vn dx ’ 1 = (1 ’ N p )5 n = 

(l,Np), and Np is the number of pressure nodes in an element as before. Solving 

for pressure, p, in equation (3.90) and substituting the result into equation (3.86) 
yields penalized discrete Navier- Stokes equation to be solved for nodal velocities 
only. The penalized discrete Navier- Stokes equations are given by 
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(3.91) 


Equation (3.91) is (2N by 2N ) system of equations, and the necessity for solving 

the conservation of mass equation and the pressure variable has been eliminated. If 
desired, the pressure field can be recovered by solving equation (3.90) in the post 
process. One advantage of the penalty method is that the number of unknowns are 
reduced so that a slight computational efficiency is achieved. But Zienkiewicz [110] 
has pointed out that solving the standard discrete Navier- Stokes equations can achieve 
comparable computational efficiency if the frontal solution technique [111] is used and 
the pressure variables are eliminated at element levels. But in this case, the pressure 
needs to be discontinuous across element interfaces. The standard solution method 
would be preferred in the flow- solid interaction analysis since the pressure solution is 
required for most of the cases. A Petrov- Galerkin finite element solution of two- 
dimensional steady laminar flow through a tube bank is shown in Figure 17, where 
the Reynolds number is 500. 


Much progress for the Petrov- Galerkin type analysis methods has been made 
recently for unsteady two-dimensional flows. Among these are: the Taylor- Galerkin 

method due to Zienkiewicz et al. [105], Donea [112], Lohner et al. [113], and Donea 
et al. [114], in which Taylor series expansion of the time derivative term is made 
prior to the spatial discretization, and the Taylor series expanded Navier-Stokes equa- 
tions are solved by the Bubnov-Galerkin method; and direct extension of the Petrov- 
Galerkin method for the two-dimensional steady flows into unsteady two-dimensional 
flows by considering the time-coordinate as an extra spatial coordinate due to Yu and 
Heinrich [105]. It is believed that most of the upwinding techniques for two- 
dimensional unsteady flows can be extended to three-dimensional unsteady flows. 
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Figure 17. A Petrov-Galerkin finite element solution of a two-dimensional 
steady laminar flow through a tube bank, (a) Pressure contour 
lines, (b) Velocity vectors, (c) Stream lines. 

Finite element computations of turbulent, separated flows also begin to appear 
in the literature. Among these are Taylor et al. [122], Smith [123], and Benim and 
Zinser [124]. In these papers, the standard k-e turbulence model for high Reynolds 
number flows has been used. In Taylor et al. [122], the flow domain was discretized 
using quadratic quadrilateral isoparametric elements as shown in Figure 18a. The 
numerically obtained mean velocity and the turbulence intensity profiles are compared 
with experimental data [125] as well as the finite difference solutions [125] in Figures 
18b and 18c. It can be seen that both the finite element method and the finite 
difference method predict the mean velocity profile accurately; however, it can also 
be seen that the finite element method predicts the turbulence intensity profile better 
than the finite difference method for this particular case. 

The boundary element method, which has been tested for incompressible laminar 
flows [127], also exists. In this method, the residual function, defined in the same 
way as in the finite element methods, is weighted against the test function and inte- 
grated over the domain. But in this method, the test fucntion is the particular solu- 
tion of the governing differential equation, and hence the domain integration is reduced 
into a boundary integration for Laplace's equation. But still one needs to integrate 
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<c) 

o: EXPERIMENTAL DATA [125] , x: FINITE ELEMENT COMPUTATIONAL RESULTS - — ■ 
FINITE DIFFERENCE SOLUTION 


Figure 18. Flow over a backward-facing step, (a) Finite element 
discretization, (b) Velocity profile, (c) Turbulent 
kinetic energy profile. 
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over the domain for Poisson type equations. Originally, the method was applicable 
only when the particular solution for the governing differential equation is available. 
For nonlinear problems, the nonlinear terms can be considered as forcing functions in 
the iterative solution procedure. The method is particularly suitable for problems 
exhibiting singularities as those associated with fracture mechanics. But the use of 
the method for turbulent flows is not certain, not to mention the flow- solid 
interactions. 


3.5.3 Body-Fitted Grid Finite-Difference Analysis of Steady Flows Through Tube Banks 

In the body-fitted grid finite -difference method, the flow domain is mapped onto 
an orthogonal computational grid; the flow equations on the computational grid are 
obtained by transforming the conservation of mass equation, the Navier- Stokes equa- 
tions , and the energy equation using the coordinate transformation relationship . The 
resulting system of partial differential equations are solved on the computational grid 
system . 

If the physical flow domain is complicated or highly curved, then the numerical 
results may not be accurate due to large truncation error involved in the grid trans- 
formation relationship. This difficulty was overcome in Kwak and Chang [106] and 
Chang and Kwak [107] by dividing the flow domain into several sub-domains which 
overlap at the interface regions. Thus, grid transformation error can be reduced and 
slight computational efficiency can be achieved by reducing the size of the system of 
equations to be solved. The body- fitted grid finite- difference method proved to be 
quite successful for laminar flows, but its extension to turbulent flows is not so 
obvious as yet . It may be expected that the truncation error in grid transformation 
may accumulate in an iterative solution procedure of coupled flow and turbulence equa- 
tions so that the computational results may not be so accurate. Success of the body- 
fitted grid finite- difference method for complicated turbulent flows may depend on 
overcoming the truncation error, one method for doing so is the sub- domain method 
[106,107]. 

In order to achieve further computational efficiency, a pseudocompressibility 
method was used in Kwak and Chang [106], Chang and Kwak [107], and Rogers et al. 

[ 108] , so that a time dependent solution method can be used to obtain a steady state 
solution. It was shown in References 106, 107, and 108 that the pseudocompressibility 
method takes less computational time than solving the steady state Navier-Stokes 
equations directly. Numerically solved grids and flow fields through a cylinder and 
two rows of cylinders, using the above described method, are shown in Figure 19 [108]. 

The body-fitted grid finite -difference method is also intended to be further 
developed for flow-solid interaction analysis of tube banks [106,107,108]. 


IV. CONCLUSIONS AND RECOMMENDATIONS 


This report provides a review of various methods used for design of tube banks 
as well as computational methods to analyze flow- solid interactions of tube banks sub- 
jected to cross flows. 

Some conclusions can be made regarding the present state of understanding of 
the flow field in tube banks and the flow- solid interactions in such tube banks. These 
are listed below. 
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Figure 19. INS3D computation of flow through cylinders, (a) Grid for a single 
post on a flat plate, (b) Grid spacing near post-plate junction for a single 
post, (c) Velocity vectors for multiple post in z/d = 0.01 plane. 

(d) Grid for calculating flow around two rows of cylinders. 

(e) Velocity vectors for single post in z/D = 0.01 plane. 



1) No major experiment has ever been made for the flow field in the region 
between tube banks and all the understanding of the flow field depends on conjectures 
by researchers in this area. 

2) Analytical models for flow -solid interaction in tube banks are based on too 
many restricting assumptions so that these analytical models cannot fully describe the 
practical problems of interest. Consequently, practical design of tube banks mostly 
depend on experimental data. 

3) Most of the analytical models use the Euler beam equation to represent thin 
cylinders. The validity of modeling thin cylinders by the Euler beam equation need 
to be studied further if these models are going to be extended to include large 
deformations and material nonlinearity. 

4) The mathematical model developed by Chen [45] is one of the most complete 
models relating to this subject , in the sense that it includes multiple number of 
vibration inducing mechanisms. Many of the coefficients, such as the fluidelastic 
stiffness coefficients used in the model, need to be supplied by experiments or by 
computational fluid dynamics as suggested by Chen [45]. But the model can predict 
only weakly coupled flow- solid interactions. Its usefulness for predicting failure 
mechanisms such as cyclic straining or hysteretic effects is limited in its present 
form because of the reasons stated in (3). 

5) Detailed numerical modeling of the flow field in the region between tube 
banks has just begun and the numerical solutions obtained until now are for steady, 
constant-viscosity flows only. 

Solving flow-solid interaction problems requires, in general, capabilities to 
solve fluid flow problems, structural dynamics problems, and coupling the flow- 
structure interface conditions. Therefore, the state of development for turbulence 
models, numerical methods such as the finite element method and the finite difference 
method, nonlinear structural analysis methods, and existing analytical methods applic- 
able to laminar as well as to turbulent flow- solid interactions have been included in 
the discussions in order to evaluate the possibility of initiating a development of 
major numerical analysis method for cross flows. 

It was found that the existing numerical analysis methods for flow-solid inter- 
actions are designed for specific problems of each researcher's interest. At this 
time, it is not known whether these methods can be naturally extended to flow-solid 
interactions for a nest of cylinders subjected to turbulent flow at high Reynolds 
numbers. 

The development stage , advantages , and disadvantages of the various numerical 
methods discussed in this report are summarized in Table 4 as they are related to the 
flow-solid interaction analysis. 

Based on these observations, some recommendations can be made. These are 
listed below. 

1) Presently, most of the structural analyses have been based on the finite 
element method. Recent developments in the Petrov-Galerkin finite element method 
enables the method to be applciable to high Reynolds number flows which in the past 
had been solved mostly by the finite difference method. It is, therefore, recommended 
that the use of the finite element method for flow- solid interaction analysis in tube 
banks be pursued further. 
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TABLE 4. COMPARISON OF DIFFERENT NUMERICAL METHODS FOR THE 
NUMERICAL ANALYSIS OF FLOW-SOLID INTERACTIONS 

(Smaller number denotes higher rating) 

CV-FDM 1 BFG-FDM 2 CV-FEM 3 FEM 4 BEM 5 FDM 6 

Capability to handle 2 3 2 1 4 1 

nonlinear equations 

Capability to handle 3 3 3 2 1 3 

complicated geometry 

Number of grids 
required 

Experience gained 
for turbulent flow 
analysis 

Experience gained NI NI NI 12 2 

for structural 
analysis 

Experience gained NI NI NI 1 NI 2 

for flow- solid 
interaction 

Computer time to 1 1 2 3NI1 

solve the same 
turbulent flows 

Computer time to NI NI NI 12 2 

solve the same 
nonlinear structural 
problems 

1. CV-FDM — Control-Volume Based Finite- Difference Method 

2. BFG-FDM — Body-Fitted Grid Finite Difference Method 

3. CV-FEM — Control- Volume Based Finite Element Method 

4. FEM — Finite Element Method 

5. BEM — Boundary Element Method 

6. FDM — Finite Difference Model 

7. NI — No Information Available. 


3 3 3 2 1 3 

1 2 3 2 NI 7 1 
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2) One of the best performing elements for shell structures are the Ahmad 
element [55] and the semiloop element [126]. Use of these elements for geometrically 
and/or materially nonlinear structures has been reported elsewhere, including Refer- 
ences 55 and 126. But use of these elements for vibration problems which require 
rigorous mathematical modeling of elasto- viscoplasticity [93], such as flow-solid inter- 
action analysis of tube banks, has not appeared as yet. It is recommended that these 
two types of shell elements need to be further tested and developed for the present 
problem . 

3) Many of the structural elements including both of the Ahmad element and 
the semiloop element are based on the eight-node serendipity element. Therefore, use 
of the serendipity element for flow calculations is recommended. The serendipity 
element or the complete quadratic element can approximate the curved and smooth 
boundary far better than linear elements. 

4) Flow-solid interaction in a nest of thin cylinders subjected to cross-flows 
can be characterized by the large deformation of the structures which entails con- 
tinuously changing flow domain. It would be advantageous, therefore, to use the 
Lagrangian formulation for the structural analysis and to use the Euler -Lagrangian 
formulation, introduced in Section 3.4, for the turbulent flow. 

5) The body-fitted grid finite -difference method has already made significant 
progress toward flow-solid interaction analysis in tube banks. It is recommended that 
this method be further developed because it can be used together with the finite 
element structural analysis codes. The only drawback envisioned with this approach 
would be the complexities associated with a computer program structure originating 
from combining two numerical methods which are quite different in their nature. 

6) Regardless of whether the finite element or the finite difference method is 
used for flow- solid interaction analysis in tube banks, difficulties may be encountered 
as a result of the use of unsteady turbulence models which are not well established 
as yet. There exist a few turbulence models well established for steady state turbu- 
lence models as has been discussed in Section 3.3. There also appeared a multiple 
scale k- e turbulence model proposed by Hanjalic et al. [120] and an extended multiple 
scale turbulence model by Kim [12]. At this time, it is not clear which 
turbulence model would be most suitable to be extended to unsteady turbulent flows. 

It should be noted that the straight forward use of the standard k- e turbulence model 
to unsteady flows has almost always been unsuccessful. Accordingly, any turbulence 
model to be used for flow- solid interaction needs to be tested against available experi- 
mental data. It is also recommended that tests for unsteady, turbulent flows be made 
of as many turbulence models as possible. 

A successful numerical modeling of the flow- solid interaction would enhance the 
efficiency and confidence in the design of tube banks. For example, fretting of 
cylinders due to large deformation can be predicted directly , viscoplastic deformations 
which will shorten the lifetime of tube banks can be predicted directly , and any 
instability of tube banks subjected to cross- flows could be predicted using the non- 
linear stability criterion analysis method [55]. 
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